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Abstract. We present a computer assisted proof of existence of a family of Lyapunov orbits 
which stretches from L2 up to half the distance to the smaller primary in the Jupiter-Sun pla- 
nar restricted circular three body problem. We then focus on a small family of Lyapunov orbits 
with energies close to comet Oterma and show that their associated invariant manifolds intersect 
' transversally. Our computer assisted proof provides explicit bounds on the location and on the angle 

C^l of intersection. 
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Xjy • 1. Introduction. The Planar Restricted Circular Three Body Problem (PRC- 

3BP) has been extensively studied throughout literature. The model has applications 
in space mission design [11, 12], explains symbolic dynamics phenomena observed in 
trajectories of comets [15] and can be used for study of diffusion estimates [13, 14]. 
All the above mentioned are associated with dynamics along invariant manifolds of 
the system. In this paper we discuss how existence of such manifolds can be proved 
within explicit bounds using rigorous-computer-assisted techniques. 

We focus on dynamics associated with the fixed point L2, its associated center 
manifold and stable/unstable manifolds. The problem has been studied by Llibre, 
^\ ' Martinez and Simo [16] where under appropriate conditions on parameters of the 

C*") . system existence and intersections of such manifolds has been proved analytically. 

In the work of Koon, Lo, Marsden and Ross [15] such invariant manifolds and their 
associated symbolic dynamics have been used to numerically explain a peculiar tra- 
jectory of the comet Oterma in the vicinity of Jupiter. Such symbolic dynamics 
has later been proved using rigorous-computer-assisted computations by Wilczak and 
Zgliczyhski [19, 20]. The work presented in this paper can be viewed as an extension 
of last-mentioned. Results [19, 20] were obtained using purely topological arguments. 
They focus on homoclinic and heteroclinic tangle between periodic orbits, without 
the detection of the manifolds themselves or angles of their intersections. Here we 
address these issues. 

In this paper wc shall first present a method for detecting of families of Lyapunov 
orbits in the PRC3BP. It is designed as a tool for rigorous-computer-assisted proofs. 
Wc apply the method to obtain a family that spans up to half a distance between the 
fixed point Li and the smaller primary in the Jupiter-Sun system. This is our first 
main result, which is stated in Theorem 3.2. The method is based on a combination 
of interval Newton method and implicit function theorem. 

We then consider a small family of Lyapunov orbits with energies close to the 
energy of comet Oterma. We prove that the family is normally hyperbolic, and give a 
tool for obtaining rigorous bounds for its unstable and stable fibers. The tool is based 
on a topological approach combined with a parameterization method. We then show 
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how fibers can be propagated to prove transversal intersections between stable and 
unstable manifolds of Lyapunov orbits. We investigate an intersection associated with 
manifolds which span from the Lyapunov orbit and circle around the larger primary. 
We obtain explicit bounds on the location of intersection and also on its angle. This 
is the second main result of the paper, which is stated in Theorem 4.1. 

Both methods which we propose are tailor made for the PRC3BP. We make use of 
the preservation of energy and reversibility of the system. Thanks to this our rigorous 
bounds for the investigated manifolds are quite sharp. 

For our method we also develop a more general tool which can be applied for 
the detection of unstable/stable manifolds of saddle - center fixed points. It is a 
generalization of the wok of Zgliczyhski [21]. This is the subject of section 6. 

The paper is organized as follows. Section 2 includes preliminaries which give 
an introduction to the PRC3BP, the interval Newton method, and introduce some 
notations. In section 3 we present a method for detection of families of Lyapunov 
orbits and apply it to the Jupiter-Sun system. In section 4 we outline the results 
for the intersections of invariant manifolds, which are then proved throughout the 
remainder of the paper. In section 5 we show how to prove that Lyapunov orbits 
are hyperbolic and foliated by energy. In section 6 we give a topological tool for 
detection of unstable manifolds of saddle-center fixed points. The method is then 
combined with parametrization method in section 7 to obtain rigorous bounds on the 
intersections of invariant manifolds. Sections 8, 9 and 10 contain respectively closing 
remarks, acknowledgements and the appendix. 

2. Preliminaries. 

2.1. The Planar Restricted Circular Three Body Problem. In the Planar 
restricted circular three body problem (PRC3BP) we consider the motion of a small 
masslcss particle under the gravitational pull of two larger bodies (which we shall 
refer to as primaries) of mass [i and 1 — fx. The primaries move around the origin on 
circular orbits of period 2ir on the same plane as the massless body. In this paper we 
shall consider the mass parameter /i = 0.0009537, which corresponds to the rescaled 
mass of Jupiter in the Jupiter-Sun system. 

The Hamiltonian of the problem is given by [1] 

H(q,p,t) = — 



2 n(t) r 2 (t)' 

where (p, q) = (qi, q2,Pi, P 2) are the coordinates of the massless particle and r\(t) and 
T2 (t) are the distances from the masses 1 — fx and [i respectively. 
After introducing a new coordinates system {x,y,p x , Py ) 

x = qi cost + q2 sinf, p x = p\ cost + P2 sin t, (2 1) 

y = — qi smt + q 2 cost, p y = — pi sin t + P2 cos t, 

which rotates together with the primaries, the primaries become motionless (see Fig- 
ure 2.1) and one obtains [1] an autonomous Hamiltonian 

H(x,y, Px , Py ) = (P* + y) 2 + (Py~*) 2 _ n(x ^ (2 . 2) 
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where 

x 2 + y 2 1 — n [i 



Cl(x,y) 



n r 2 



n = \J (x - /i) 2 + y 2 , r 2 = \/(x + 1 - fi) 2 + y 2 . 
The motion of the particle is given by 

q = JVH(q), (2.3) 

where q = (x, y,p x ,Py) € K 4 , J = ( ^ ^0 ) ^ * S & dimensional identity 
matrix. 

The movement of the flow (2.3) is restricted to the hypersurfaces determined by 
the energy level h, 

M{h) = {(x,y,p x ,p y ) g M. 4 \H{x,y,p x ,p y ) = h}. (2.4) 

This means that movement in the x, y coordinates is restricted to the so called Hill's 
region defined by 

R(h) = {(x,y)€R 2 \n(x,y) > -h}. 




The problem has three equilibrium points L\^Li^L^ on the x-axes (see Figure 
2.1). We shall be interested in the dynamics associated with L 2 , and with orbits of 
energies higher than that of L 2 . The linearized vector field at the point L 2 has two 
real and two purely imaginary eigenvalues, thus by the Lyapunov theorem (see for 
example [16]) for energies h larger and sufficiently close to H(L 2 ) there exists a family 
of periodic orbits parameterized by energy emanating from the equilibrium point L 2 . 
Numerical evidence shows that this family extends up to and even beyond the smaller 
primary /j, [3]. 
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Fig. 2.2. Considered by us family of Lyapunov orbits in green (spanning between two orbits in 
blue), together with the Lyapunov orbit for the energy of the comet Oterma h — 1. 515 in red. 

The PRC3BP admits the following reversing symmetry 

S(x,y,p x ,p y ) = (x, -y, -p x ,p y ). 

For the flow <p(t, q) of (2.3) we have 

S(<f>(t t q)) = <t>(-t,S(q)). (2.5) 

We will say that an orbit q(t) is 5-symmetric when 

S(q(t))=q(-t). (2.6) 

Each Lyapunov orbit is S-symmetric. It possesses a two dimensional stable man- 
ifold and a two dimensional unstable manifold. These manifolds lie on the same 
energy level as the orbit and their intersection, when restricted to the three dimen- 
sional constant energy manifold (2.4), is transversal. These invariant manifolds are 
S'-symmetric with respect to each other, meaning that the stable manifold is an image 
by S of the unstable manifold (see Figure 2.3 for the unstable manifold, and Figure 
2.4 for the intersection of manifolds). All these facts are well known and extensively 
studied numerically. 

Our aim in this paper will be firstly to provide a rigorous-computcr-assisted proof 
of existence of the manifold of Lyapunov orbits over a large radius from L 2 (see Figure 
2.2). Secondly, using rigorous-computer-assistcd computations, we shall show that for 
orbits with energies close to the energy of comet Oterma h = 1. 515 their associated 
stable and unstable manifolds intersect transversally. Even though such intersections 
are well known from numerical investigation, to the best of our knowledge this is a 
first rigorous proof of their existence. 

2.2. Interval Newton Method. Let X be a subset of E™. We shall denote by 
[X] an interval enclosure of the set X, that is, a set 

[X}=ir? =1 [a i ,b i ]cM n , 

such that 




x c [X]. 
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Fig. 2.3. The Lyapunov orbit in red, its unstable manifold in green, and the intersection of the 
unstable manifold with section {y = 0} in blue, projected onto x,y,p x coordinates. The figure is for 
the energy of comet Oterma h = 1.515 in the Jupiter-Sun system. 
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Fig. 2.4. The Lyapunov orbit in red, its unstable manifold in green, stable manifold in purple, 
and their intersections with section {y = 0} in blue, projected onto x,y coordinates (left) and x,p x 
coordinates (right). The figure is for the energy of comet Oterma h = 1.515 in the Jupiter-Sun 
system. 



Let / : K" -> K" be a C 1 function and U C K". We shall denote by [Df(U)} the 
interval enclosure of a Jacobian matrix on the set U. This means that [Df(U)] is an 
interval matrix defined as 



[Df{U)\ 



AG 



inf — (a;), sup — (x) 
xeu axj x ^u axj 



for all i, j = 1, . 
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Theorem 2.1. [2] (Interval Newton method) Let f : M" — >• W l be a C 1 function 
and X — n™ =1 [aj,6 i ] with < hi. If [Df(X)\ is invertible and there exists an xq in 
X such that 

N(x , X) := xo - [Df(X)}- 1 f(x ) C X, 

then there exists a unique point x* £ X such that f(x*) = 0. 

2.3. Notations. Throughout the paper we shall use a notation 4>(t,x) for the 
flow, and $t(x) = <j>(T-> %) for a time T shift along trajectory map of (2.3). For points 
p = (x, y) we shall write n x p and ir y p, to denote projections onto coordinates x and y 
respectively. We shall also use the following notation for a cartesian product of sets 
n™ =1 C/i = U\ x . . . x U n . For A, B C K" we shall use a notation A + B = {a + b\a £ 
A,b£ B}. 

3. Existence of a Family of Lyapunov Orbits. In this section we shall 
present a method for proving existence of Lyapunov orbits far away from L^. The 
result is in the spirit of the method applied by Wilczak and Zgliczyhski in [19, 20] 
for a Lyapunov orbit with energy h = 1. 515 of the comet Oterma. Our result differs 
from [19, 20] by the fact that we obtain a smooth family of orbits over a large set, 
whereas in [19, 20] a single orbit was proved. 

We shall consider orbits starting from points of the form 0, 0,p y ) with x inside 
an interval 



Ix \ I'X ; Ix\ • 

w [-0.96602315,-0.933] C 



+ At - 0.933), -0.933 



(3.1) 



Since 7r T L 2 ~ —0.93237 we see that 4 < i(-l + /j — n x L2), so the interval I x 
stretches from half the distance between the smaller primary and L 2 , almost up to 
(see Figure 2.2, where the orbits are depicted in green, and stretch between an inner 
and outer orbit depicted in blue). 

Let us consider a section E = {y = 0} and a Poincare map P : S — > S of (2.3). 
We shall interpret the Poincare map as a function from M 3 to M 3 with coordinates 
x,p x ,p y . If for a point q = (x, 0,p y ) £ E we have 7r Px P(g) = 0, then by the symmetry 
property (2.5) the point q lies on a periodic orbit (the Poincare map P makes a half 
turn along the orbit starting from q). 

Let us introduce the following notation 

/ : R 2 E, 



f(x,p y ) = 7T Px P(a;,0,Py). 
To find a periodic orbit for some fixed x it is sufficient to find a zero of a function 

9x(p v ) ■= f(x,p y ). 

Let DP = {dPij) i j_ 1 2 3 be the derivative of the map P, with indexes 1, 2,3 corre- 
sponding to coordinates x,p x ,p y respectively. 

Lemma 3.1. Let I and Ji for i = 0, 1 be closed intervals such that Jo, J\ have 
the same center point p® and Jq C J\. Let x° be the center point of I. Let agl and 
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Fig. 3.1. The bound for a curve of points q(x) = (x,0,0, k(x)) on Lyapunov orbits. 



U , U C £ 



be sets defined as (see Figure 3.1) 
U = {x°} x {0} x J , 

U = {(x,0,p v )\x G I,p y = a {x — X ) + I, I G Ji) 



7r Px P(x° ; 0,pg) 
dP(i7o)2 3 



C Jo, 



and 



a\ < U| (| Ji| - | Jq|) /oraIIae[a,a]:= 



dP(C7) 



2 1 



(3.2) 



(3.3) 



(3.4) 



dP(U) 23 

then there exists a smooth function k : I — > R smc/i i/iai /or any a; G / a point q(x) = 
(x, 0, 0, k(x)) lies on an S- symmetric periodic orbit of (2.3). Moreover, k'(x) G [a, a] 
and q(x) G U for all x G /. 

Proof. Existence of a unique point k(xq) G Jo for which g Xo (K,(xo)) = follows 
from (3.3), which implies 

pi - [DgMT 1 g xo (p° y ) C N C J , 

combined with interval Newton method (Theorem 2.1). 

For (3.4) to hold we need to have ^ dP(U)23- For (x,0,p y ) G U we have 
J^-(x,p y ) G dP(U)23 hence -^-(x,p y ) ^ 0. This means that we can apply the 
implicit function theorem to obtain a curve k(x) for which f(x,n(x)) = 0. We now 
need to make sure that the curve k is defined on the entire interval /. At each point x 
for which (x, 0, n(x)) G U is defined, by the implicit function theorem we know that 

%(x,k(x)) 



k'{x) 



dP(U) 



2 1 



dP(U)i 



This, by assumption (3.4), means that we can continue the curve from k(x°) to the 
whole interval / (see Figure 3.1). □ 

To apply Lemma 3.1 we first compute numerically a sequence of points (see Figure 

3.2) 

^ = (^,0,0,^,) for * = (),..., 15 000, 



8 



MACIEJ J. CAPINSKI 



-0.76 
-0.78 
-0.8 
-0.82 
-0.84 

> 

Q. 

-0.86 
-0.88 
-0.9 
-0.92 
-0.94 

-1 -0.98 -0.96 -0.94 -0.92 

x 

Fig. 3.2. Numerical plot of k(x), consisting of 15 000 points q® on Lyapunov orbits (in red). 
The point L2 is in green. The blue line x = —1 + fi gives an indication of the position of the smaller 
primary along the x coordinate. 

where I x , I x are defined in (3.1). The q® are non-rigorously, numerically computed 
points on Lyapunov orbits. We then compute (non-rigorously) a sequence of slopes 
(see Figure 3.3) 

Oj 68 i = 0,..., 15000, 

define 

r = -J— -(4 -I x ) w 10~ 6 • 1.1007716, 
15 000 2 v — 

h = x°i + [-r, r) , 
Jo,i=P^ + 10- 13 -[-l,l], 
J u =pO i + 10- 8 - [-5,5], 

and consider sets 

U = {x a t } x {0} x J 0ii , 

Ui = {(x,0,p y )\x G Ii,p y =ai(x— x°) + i,i G Ji ti } . 

We apply Lemma 3.1 repeatedly 15 000 times, and obtain the following theorem. 

Theorem 3.2 (First main result). Let I x be the interval from (3.1). Then there 
exists a curve q(x) = (x, 0, 0, k(x)) of points on Lyapunov orbits with k : I x — > R, 
which lies within a 5 ■ 10 -8 distance from the piecewise linear curve joining the 15 000 
points q® on Figure 3.2. 

The proof of Theorem 3.2 took 5 hours and 43 minutes on a standard laptop. 

Remark 3.3. Using above described method it is impossible to continue with the 
orbits to L2- At the fixed point one would need to apply alternative methods, such as 
the method of majorants [18], Lyapunov theorem by tracing the radius of convergence 
of the normal form [l 7], or topological- computer- assisted tools such as [5, 7j. 

4. Outline of Results for Intersections of Invariant Manifolds. In the 

reminder of the paper we shall focus our attention on orbits starting from q(x) = 
(x, 0, 0, k(x)) with x G / for 

I = [LJ] ■= x° + [-1,1] • 10~ 9 , (4.1) 
x° = -0.9510055339445208. 
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Fig. 3.3. Numerical plot of k'(x), consisting of 15 000 points a; 

Such orbits have energy close to the energy of the comet Oterma h = 1.515. 

Let us introduce a notation A for a family of Lyapunov orbits, which start from 
q(x) with x G I 

A = {4>(t,q(x))\t G R,q(x) = (x, 0, 0, k(x)), x G 7}. (4.2) 

For x € I, let L(x) C A denote the Lyapunov orbit which starts from q(x). 

Throughout the reminder of the paper we shall prove the following theorem. 

Theorem 4.1 (Second main result). A is a normally hyperbolic invariant man- 
ifold with a boundary. Each orbit L(x) C A possesses a two dimensional stable man- 
ifold W s {L(x)) and a two dimensional unstable manifold W u (L(x)). The manifolds 
W 8 (L(x)) and W u (L(x)) intersect and the intersection, when restricted to the con- 
stant energy manifold M(H(L(x))), is transversal (see (2.4) for definition of M). 

Numerical plots of the intersection of manifolds that we shall prove are given in 
Figure 2.4. 

Theorem 4.1 will be proved with computer assistance. During the proof we shall 
obtain rigorous bounds on the region and the angle at which the manifolds intersect 
(see Figure 7.6). 

The size of interval I (4.1) is very small. When translated the real life distance 
in the Jupiter-Sun system, its length is just slightly over one and a half kilometer. 
This is practically a single point. We need to start with such a small set to obtain 
our result. Thanks to this we obtain sharp estimates on the intersection of W s (L(x)), 
W u (L(x)). To consider a larger set of Lyapunov orbits one would need to iterate the 
procedure a number of times. This can be done without any difficulty apart from 
necessary time for computation. The proof of Theorem 4.1 took 46 minutes on a 
standard laptop. Using clusters one could cover a larger interval I in reasonable time. 

5. Hyperbolicity of Lyapunov Orbits and Foliation by Energy. In this 
section we shall show that each orbit L(x) C A lies on a different energy level. We 
shall also show that each orbit L(x) (when considered on its constant energy manifold) 
is hyperbolic. In other words, we shall show that A is a normally hyperbolic manifold 
with a boundary. 

We start with a simple remark. 

Remark 5.1. If for all x G I we have -4-H(q(x)) 7^ 0, then Lyapunov orbits 
with different x have different energies. Note that the set U and the bound on the 
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derivative of n'(x) from Lemma 3.1 can be used to obtain 



We shall now give a simple lemma which can be used to show that our Lyapunov 
orbits arc hyperbolic. 

In what follows in this section, let P : X — > £ be a second return Poincare map 
for £ = {y = 0}. This means that each point q{x) = (x, 0, 0, k(x)), with x £ I, is a 
fixed point of P. We shall interpret the Poincare map as a function from M 3 to K 3 
with coordinates x,p x ,p y . 

Lemma 5.2. Let U be the set given by (3.2) in Lemma 3.1. Assume that for any 
1x2 matrix A 



Ae 



-W 

dx ) V 



dH 



dH 

dp y 



(CO 



and any 2x2 matrix B 

B e 



dP 2 i 
dP 3 i 



.4 



dP 22 rfP 23 
dP 32 rfP 33 



(CO 



(5.1) 



(5.2) 



t/ie spectrum of B consists of two real eigenvalues Ai,A2 satisfying |Ai| > 1 > | A2 1 . 
Then for any x G / the Lyapunov orbit starting from q(x), restricted to the constant 
energy manifold M(H(q(x))), is a hyperbolic orbit. 

Proof. Let us fix some x € I. For our assumptions to hold, A from (5.1) needs to 



be properly defined. This means that ^-{q(x)) 7^ 0. By the implicit function theorem 



there exists a function x{p Xl p y ) with x(0, n(x)) 
H(q(x)) and 



x such that H(x(p x ,p y ),0,p x ,p y ) 



dx dx 
dp* dp y 





dH. 


OH \\ 




dpx 


dPy J j 


\ dx 







(5.3) 



The Lyapunov orbit starting from q{x) is contained in the constant energy man- 
ifold M(H(q(x))). Let us consider V = M(H(q(x))) (1 {y = 0} and a Poincare map 
P : V — > V. In a neighborhood of q(x) the manifold V can be parameterized by 
(p x ,Py) ■ Since 



P(Px,Py) = ■K(p x ,p v )P(x(Px,P V ),Px,Pv) 



we have 



DP(0,k(x)) 



(5.4) 



9P \ f Qx_ Qx_ 

7T (Pa 1 ,P v ) Q x J \ dp* dp y 



dP 22 dP 23 
dP 32 dP 33 



(x,0,k(x)). 



By (5.3), (5.4) and our assumption about the spectrum of B of from (5.2), follows that 
(0,k(x)) is a hyperbolic fixed point for the map P. This means that the Lyapunov 
orbit starting from q(x), restricted to the constant energy manifold M(H(q(x))) is 
hyperbolic. □ 

Remark 5.3. Since B from (5.2) is a 2 x 2 matrix, estimation of its eigenval- 
ues is straightforward. Here we profit from the the reduction of dimension made by 
restricting to a constant energy manifold. 
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Since we consider a small part of the family of orbits (4.1), we can obtain a much 
tighter enclosure of the curve k(x) for x G I than from Theorem 3.2. Let 



p o = -0.836804179646973 J = p° y + [-1, 1] • 10 

v 

and 



4.506866203376769 Ji=p° v + [-1, 1] • 10" 12 



(5.5) 



U = {(x,0,0,p y )\x G I,p y = a(x-x°) +i.,ie Ji} . (5.6) 

PROPOSITION 5.4. For x G J, wit/i I from (4-1), we have q(x) = (x,0, 0,k(x)) C 
[/ and 

k'(x) g [-4.506980818,-4.506751634], (5.7) 

-^-H(q(x)) G [-0.3670937615,-0.3670674516], (5.8) 
ax 



if(I,0,0,a(I-x o ) + Ji) G [-1.514999999635,-1.514999999631], 
#(7,0,0,a(7-x ) + Ji) G [-1.515000000369,-1.515000000365]. 

Moreover, the orbits (when considered on their constant energy manifolds) are hyper- 
bolic, and we have following bounds for the eigenvalues 

Ai G [1450.24, 1481.68], (5.9) 
A 2 G 10~ 4 [6.74909,6.89541]. 

Proof. The proof was performed with computer assistance. It required no subdi- 
vision of U and the computation took less than two seconds on a standard laptop. 

Existence of q(x) C U was shown using Lemma 3.1. From it also follows the 
bound (5.7) for k'(x). The bound (5.8) follows from Remark 5.1. Hyperbolicity and 
bounds (5.9) follow from Lemma 5.2. □ 

6. Cone Conditions and Bounds for Unstable Manifolds of Saddle- 
Center Fixed Points. In this section we provide a topological tool that can be 
used for rigorous-computer-assisted detection of unstable manifolds of saddle-center 
fixed points. The method is a modification of [21], where instead of saddle-center 
a standard hyperbolic fixed point was considered. The content of this section is a 
general result. In section 7 we return to the PRC3BP and show how to apply it for 
the proof of Theorem 4.1. 

Let F : W n — s- M" be a C k diffeomorphism with a fixed point v* G M ra and k > 1. 
Assume that for eigenvalues Ai, A2, . . . , A„ from the spectrum of DF(v*) we have 

|rcAi|>m>l, (6.1) 
|reAi I < m for i = 2, . . . , n. 

Let W u (v*) denote the unstable manifold of v* associated with the eigenvalue Ai 

W u (v*) = {v\ \\F- n (v) ~v*\\< Cm-" for all n G N and some C > 0} . 

Let u = 1 and c = n — 1. The notations u and c will stand for "unstable" and 
"central" coordinates of F at v* . Consider two balls B u and B c , of dimensions u and 
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c respectively, such that B u x B c is centered at v* . For a point iiGlR" x K c we shall 
write v = (x, y) , with x G M. u , yGK c . In these notations we shall also write the fixed 
point as v* = (x* , y* ) . 

Remark 6.1. We do not need to assume that (x, 0) is the eigenvector associated 
with Ai and that vectors (0,y) span the eigenspace of A 2 , . . . , A„. For our method to 
work it is enough if these vectors are "roughly" aligned with the eigenspaces. This is 
important for us, since in any computer assisted computation it is usually not possible 
to compute the eigenvectors with full precision. 

Let a G K, a > and consider a function Q : R u x M c -> K 

Q(x,y) = ax 2 - ||y|| 2 . 

For vo G K u x R c we shall use a notation Q + (vo) for a cone 

Q+(t>o) = {u|Q(i;-vo)>0}. 

Let us assume that a is chosen sufficiently small so that Q + (v*) n B« x B c does not 
intersect with B u x <9£? c (See Figure 6.1). 

Definition 6.2. We shall say that h : B u — > B u x _B C is a horizontal disc in 
B u x B c for cones given by Q if h(x*) = v* , 7r x /i(x) = x and for any xi 7^ X2 holds 
Q (/i(xi) - ft(x 2 )) > 0. 



«0\ 












/ 







/12 
'/l 



X() 



(x,u; u (x)) x 

Fig. 6.1. Construction of the curve (x, wj"(x)) which lies on the unstable manifold of v* . 

Lemma 6.3. Assume that for any v\,V2 G suc/i i/iai Q(i>i — V 2 ) > 0, we 

have 

Q(F( Vl ) - F(v 2 )) > 0. (6.2) 

Lei m &e </ie constant from (6.1). If for any uGfi, »/ u*, Q(w — t>*) > holds 

\\F{v)-v*\\ >m\\v-v*\\, (6.3) 

i/ien W u (v*) C Q + (u*). Moreover, there exists a function w u : B u — > £? c such that 
(id,w u )(B u ) = W u (v*) n U, and for any Xi,x 2 G B u , xi 7^ x 2 

Q(( Xl , w " ( Xl ))-(x 2 , w u (x 2 )))>0 

and 

||(w u )'(x)|| < y/a forallxeBu. 



(6.4) 
(6.5) 



Proof. We shall first show that for any xo G B u \{x*} there exists a point vq = 
(xo,w u (xo)) G Q + (v*) such that w G W u (v*). Let /io(x) = (x, y*) bc a horizontal 
disc (Sec Figure 6.1). Observe that F(ho(x*)) = F(x*,y*) = v*. By assumptions 
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(6.2), (6.3) the curve F(h (x)) is contained in Q + (v*) and F{h {dB u ))C\B u x B c = 0. 
Moreover, by assumption (6.2) for any X!,x 2 € B u , xi =/= x 2 

Q(F(/i (xi)) - F(h (x 2 ))) > Q(/i (xi) - fto(x 2 )) > 0, 

which means that {F(ho(x))\x G B u } (~l B u x B c is & graph of a horizontal disc. Let 
us denote this disc by hi and observe that hi(x*) = v*. In other words, let hi be the 
graph transform of the disc ho. 

Taking F(hi(x)) and applying an identical argument, we observe that 

{F(hi(x))\xe B U }DB U x B c 

is a graph of a horizontal disc h 2 ■ Repeating this procedure we can construct a sequence 

of horizontal discs ho, hi, hi, For a fixed Xo, due to compactness of closure of B c , 

there exists a subsequence h^ixo) convergent to some point vq G B u x cl£> c . For any 
i,n G N with hi > n the point F~ n (hkt (xo)) lies on the graph of h^-n and hence is 
also in Q + (v*). This means that for any n G N 

F- n (v ) = lim F- n (h ki (x )) G Q+(w*). 

By assumption (6.3) we have 

\\F-»(vo)-v*\\ <^\\v -v*\\, 

which means that w S W u («*). By construction 7r x wo = x , hence we can define 
w u (x ) := 7r y u . 

By the stable/unstable manifold theorem, there exists a small interval I £ = (x* — 
e,x* + e) in which {(x, w u (x))\x £ I e } is a curve which gives full description 
of W u {v*). Since (x,w"(x)) C Q + (w*) we have (1, {w u )'(x*)) G Q+(0). Since for 
sufficiently small e the vector (l,(w ll )'(x)) is arbitrarily close to (1, (u>")'(x*)), for 
xi,x 2 G I e 

Q((x 1 ,w u (xi))~(x 2 ,w u (x 2 ))) >0. (6.6) 

Iterating the curve (x, w u (x)) through F, by (6.2), (6.3) we obtain our function w u : 
B u — > B c . Note that by our construction for any Xi,x 2 G B u inequality (6.6) holds. 
This implies that for any Xi , x 2 G B u 

\\w u (xi)-w u (x 2 )f 

2 < a > 

|xi - x 2 | 

which in turn gives (6.5). □ 

Remark 6.4. Lemma 6.3 can easily be generalized to higher dimension ofW u (v*). 
The proof would be identical, taking Q(x, y) = a||x|| — ||y|| 2 . Here we have set up our 
discussion so that W u {v*) is one dimensional simply because this is what we shall 
need for our application to the PRC3BP. 

Remark 6.5. By taking the inverse map, Lemma 6.3 can be used to prove exis- 
tence of stable manifolds. 

To verify assumptions (6.2) and (6.3) in practice, it is best to make use of an 
interval matrix A = [DF(Q + (v*))]. Then for any vi,v 2 G Q + {v*) we have 

F{v x ) - F{v 2 ) = f DF(v 2 +t (vi - %)) dt ■ (vi - v 2 ) G A (vi - v 2 ) . (6.7) 
Jo 
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This means that 

F{v) - v* C A (v - v*) . (6.8) 

To verify (6.3) using (6.8) we can apply Lemma 10.1 from the Appendix. 

Let us now turn to verification of (6.2). Let Cq be a diagonal matrix such that 
v t Cqv = Q(v). Equation (6.7) gives an estimate 

Q (F(«i) - F(v 2 )) C («! - v 2 f A T C Q A ( Vl - v 2 ) . (6.9) 

To verify (6.2) using (6.9) we can apply Lemma 10.2 from the Appendix. 

7. Rigorous Bounds for Invariant Manifolds associated with Lyapunov 
Orbits. In this section we give a proof of Theorem 4.1. In sections 7.1 and 7.2 
we shall show how to apply the method from section 6 to detect fibers of unstable 
manifolds of Lyapunov orbits. In section 7.3 we shall show how to prove that the 
manifolds intersect. Using these results, in section 7.4 we give a computer assisted 
proof Theorem 4.1. 

7.1. Parameterization Method. The method from section 6 requires a good 
change of coordinates which " straightens out" the unstable manifold. We shall obtain 
such a change of coordinates using a parameterization method. In this subsection we 
give an outline of this procedure. 

In this section we shall fix some x G I and show how to find an unstable fiber of 
a point 

go = q(x) = (x, 0, 0, k(x)) G L(x). 

We shall use a notation t = r(go) for the return time along the trajectory. The point 
go is a saddle center fixed point for a r-time map <E> T : R 4 — > R 4 . 

Let C denote a matrix which brings £>$ r (go) to real Jordan form. By $ T : R 4 — ► 
R we shall denote the time r map in the linearized local coordinates 

$ r (v) := C- 1 ($ T (g + Cv) - g ) . 

Let W U ($ T ,0) denote the unstable manifold of $ T at zero. If we can find a 
function 

K = (K , K 1 ,K 2 ,K 3 ) : R — s- R 4 , 

which for all x in an interval 7o = [x, x], x < < x, is a solution of a cohomology 
equation 

$ r {K(x)) = K{Xx), (7.1) 

then K(x) C W U ($ T , 0) for x G h- 

Once K is established we can consider a nonlinear change of coordinates 

i/)= (^0,^1,^2,^3) :M 4 ^M 4 

defined as 



V^o (x,yi,y 2 ,y 3 ) = K (x) - (yiK[(x) + y 2 K' 2 (x) + y 3 K' 3 (x)) , (7.2) 
ipi (x, yi, Y2, Y3) = Ki(x) + YiK' (x) for i = 1, 2, 3. 
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Note that ip(x,0) = K(x) gives points on the unstable manifold of the fixed point for 
the map $ T . The intuitive idea behind (7.2) is to orthogonalize coordinates around 
K(x) (see Figure 7.1). 

Let us define a local map 

F = V> -1 o$ T of/i. 
Such map will play the role of F from section 6. Observe that 

{C^P (K{x)) + qo\x e Io} C Cty (W u (F, 0)) + qo = W u ($ T , q ) C W u (L(x )). 



W u (F,0) 



W U ($ T ,0) 

Fig. 7.1. The nonlinear change of coordinates ip- 



7.2. Bounds for Unstable Fibers through Parameterization and Cone 
conditions. The map ip (7.2) gives us a change of coordinates which locally "straight- 
ens out" the unstable manifold. The problem with applying the procedure from sec- 
tion 7.1 in practice lies in the fact that usually finding an analytic formula for K 
satisfying (7.1) is impossible. The best that can be done is to find a K which is 
a polynomial approximation of a solution of (7.1). This can be done by expanding 
$ T into a Taylor series and inductively comparing the coefficients in (7.1) (for a de- 
tailed description of this method we refer the reader to [4]; see in particular Section 
4 and Theorem 4.1). If we find such an approximate solution of (7.1), then the set 
{(x, 0)|x <G Io} is no longer the unstable manifold for F (defined by (7.3)), but its 
approximation. Even though then our description of the unstable fiber is not entirely 
accurate, we can apply the method from Section 6 to obtain a rigorous enclosure of 
W u (F, 0). This enclosure can then be transported to the original coordinates. 

In this Section we shall assume that go is an arbitrary point close to q(x) = 
[x, 0, 0, k(x)) for x S /, C is some given matrix and K : R — s- R 4 is some given 
polynomial and that ip is defined by (7.2). 

Remark 7.1. Let us stress that the point qo is a numerical approximation ofq(x), 
the matrix C will be a (non-rigorous) numerically obtained estimate for the change 
into Jordan form of the map $ r . We do not assume that this change is rigorously 
computed. This is practically impossible due to the fact that we do not have an analytic 
formula for D$> T (qo). For us the matrix C is simply some approximation of the matrix 
which takes D$ T (qo) into Jordan form. Let us note that it is not difficult to find an 
interval matrix C^ 1 such that the inverse matrix of our C is contained in C . 

Remark 7.2. In our setting the polynomial K is an approximation of the solution 
of (7.1). In practice we cannot obtain a fully rigorous solution of (7.1). It is important 
to emphasize that we also do not have an inverse of ip. It is also not simple to find 
good rigorous estimates for the function ip -1 due to the fact that K is a high order 
polynomial. We shall therefore set up all our subsequent computations so that we will 
never need to use the inverse function ofip. 
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For x € I, let r(x) be the period of an orbit L(x) C A. We define a map 

F = i- x o^ T{x) o^. (7.3) 

Note that for each iG/wc have a different map F, We omit this in our notation for 
simplicity, and also because below methods for obtaining rigorous bounds for F and 
its derivative shall work for all x E I. 

We shall first be interested in computing rigorous bounds for F(U). It turns out 
that (7.3) is impossible to apply since we do not have a formula for ip . Even if we 
did, direct application of (7.3) in interval arithmetic would provide very bad estimates 
due to strong hyperbolicity of the map. We use a more subtle method. 

We shall first need the following notations. Let T denote an interval such that 
t(x) <G T for all x <E I. Let A £ 1 be some number close to an unstable eigenvalue 
of D<l ) T ( x )(g(x)) for some x € I. We shall slightly abuse notations and also consider 
A : R 4 — > R 4 as a function defined on v = (x, y) e R x R 3 as 

A(x,y) := (Ax,y) . 

The following Lemma allows us to obtain rigorous bounds on pre-images of F 
from (7.3). 

Lemma 7.3. Let U x c R 4 be a given set. Let G : R x M 4 x R 4 -> R 4 be defined as 
G(r, v u v 2 ) = $ T (CiP (vi) + q ) - (CiP(X(v 2 )) + q ) . (7.4) 
Let U2CR 4 be a set and A (U2) be an interval matrix defined as 

A(U 2 ) = -[CDiP(X(U 2 ))DX\. 

If 

N(T, wo, E/i, U 2 ) := v - [A i^)]" 1 [G(T, U u v )} C U 2 , (7.5) 



then 



F(Ui) C A (U 2 ) . (7.6) 

Proof. The proof is given in the Appendix in section 10.2. See also Remark 10.4 
for comments on practical application of the lemma. □ 

Remark 7.4. The choice of the function G is motivated by the following diagram. 

R 4 ^> R 4 

t C + q a t C + qo 

R 4 ^> R 4 

TO 4 



T/ie diagram is not fully commutative, hence the bracket for A. Intuitively, for v = 
(x, 0) £ R x R 3 t/ie diagram should "almost commute". Even though this statement is 
nowhere close to rigorous, it might make the method and proof of Lemma 7. 3 more 
intuitive. 
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We now turn to the computation of rigorous bounds for the derivatives of (7.3). 
For any (x, y) contained in a set B C K 4 we have the following estimates 



DF(x, y) = (£>V (^(x, y))) 1 C^D^ (Cty(x, y) + q ) CD^(x 7 y) 



(7.7) 



c 



(D^(F(B))Y 
[DF{B)\ . 



[C- 1 ] ■ [D<f T {Ci>{B) + go)] -C-lDiP (B) 



Note that to compute \DF(B)] from (7.7) we do not need to use 
Remark 7.5. Using Lemma 7.3 and (7.7) we can in practice compute rigor- 
ous bounds for [F(B)] and [DF(B)]. We perform such computations in Section 7.4-1 
with the use of CAPD library (http://capd.ii.uj.edu.pl/). The library allows for com- 
putation of rigorous estimates for $t and its derivative and for rigorous-enclosure 
operations on maps and interval matrixes. 

Proposition 5.6 gives a bound on a set U (5.6) which contains all fixed points 
q(x), with x G /, of the map & T / X y This set can be transported to local coordinates 
(x,y). Let Bq C M 4 be such set that 

{^(C-'iqix) - q ))}\x e 1} C B . 

Such set can easily be computed using for example a technical Lemma 10.5 from the 
Appendix. 

Taking a four dimensional set (see Figure 7.2) 



B 



U Q- 

v£B 



(V)C 



using (7.7) and Lemmas 10.1, 10.2 to verify assumptions of Lemma 6.3, we can obtain 
a bound for the unstable fibers of all q(x) for x £ I. The obtained bound is computed in 
local coordinates (x, y), but can easily be transported back to the original coordinates 
(x,y,p x ,p y ) of the system. Detailed results of such computation will be presented in 
section 7.4.1. 

B E 



Bl 



Fig. 7.2. Local bound on the unstable manifold. Each fixed point ip (C 1 (q(x) 
x a I , lies in Bq and its unstable manifold is contained in B = U veBo Q + (v). 



go)), for 



Remark 7.6. Let us emphasize that to apply the method it is enough to use a 
single point go, single matrix C and single nonlinear change ip. It is not necessary to 
use different changes to local coordinates for different x G I . 

Remark 7.7. Let us note that from the fact that W s (L(x)) is S -symmetric to 
W u (L(x)), without any effort we also obtain mirror bounds for fibers of W s {L(x)) . 

7.3. Transversal Intersections of Manifolds. In this section we discuss how 
the bounds for fibers of q(x) discussed in section 7.2 can be used to prove transver- 
sal intersections of manifolds W u (L(x)) and W s (L(x)) for L(x) C A (see (4.2) for 
definition of A). 
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Let x l ,x r G M. be such that x' < x r and tt x B < x l ,x r . Let B c C M 3 be such that 
ir y B C B c . Let Be, B l E , B r E be defined as (see Figure 7.2) 

B E = [x l ,x r ] x B c , 
B l E = {x 1 } x B c , 
B E = {x r } x B c> 

and let 

V + = {(x,yi,y 2 ,y 3 ) 6R 4 |x = l,y< 6 fori = 1,2, 3}, 

V- ={(x,yi,y 2 ,y 3 ) eM 4 |x = -l,y 4 e [-^3, \/a] fori = 1,2,3}, 
V = {7-yK' G F + ,7 > 0} U {7«|v eF,7> 0}. 

Note that 

Q + (0) c V. 

Consider a section 

s = { v = o} n {x > o} n {vl < 2 (H(L(x)) + n(x, y ))}. 

This shall be a section where we detect the intersection of W u (L(x)) and W s (L(x)) 
(see Figures 2.3, 2.4). Let ^ be the flow of (2.3) and define 

r(q) = inf{i > : <j)(t,q) G E}, 
: B £ -> E, 
Q (x, y) - 0(r(C7V(x, y) + 9°), C^(x, y) + q°). 

Lemma 7.8. Assume that for F defined in (7.3) assumptions of Lemma 6.3 hold. 
If also 

^G{B l E )<0, ir P *S(B E )>0, (7.8) 

then for any x £ I (with I defined in (4-1)) the manifolds W u (L(x)) and W s (L(x)) 
intersect. 

Moreover, if for any v + G V + and v~ G V~ 

ttx [DQ (B E )] v+ > 0, % Vx [DQ (B E )} v+ > 0, (7.9) 
K X [DQ (B E )} v- < 0, 7t Px [DQ (B E )} v~ < 0, 

then for each fixed x G / the intersection is transversal on the constant energy manifold 
M(H(L(x))) (see (2.4) for definition of M). 

Proof. Let us fix an x G I. First let us observe that because energy (2.2) is 
preserved, the manifold M(L(x)) n E can be parameterized by x,p x since 

Py = Py (x, Px ) = ^2(H(L(x)) + n(x, y)) -p*+x (7.10) 

is well defined. 

By Lemma 6.3 we know that in local coordinates x, y the unstable fiber of q(x) is 
a horizontal disc in B. This disc is a graph of a function w u : B u — > B c and for any 
xi , X2 € B u such that xi ^ x 2 

(x 1 ,w u (x 1 )) - (x 2 ,w u (x 2 )) G Q + (0) C V. 
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The disc also passes through the set Be (see Figure 7.2). 

In the statement of our lemma we implicitly assume that Q (x, y) is well defined 
for (x, y) € Be- This means that 

W u {L(x)) nsng (b e ) = {G (x, w u (x)) |x g [x',x r ]}. (7.11) 

Let us introduce a notation 

: [x',x'-]^M 2 , 

wl{x)=7T x , Px g(x,W U (x)). 

By (7.10) and (7.11) the curve (x) parametrizes a fragment of the intersection of 
W u {L(x)) with S. By assumption (7.8) 

TTp.^ (x») = TTp.0 (x'.^x')) G ({x<} X S c ) - (B^) < 0, 

(x r ) = ^e(x r lW "(x r )) G TTp.0 ({x'} x B c ) = tt Px Q (B r E ) > 0, 
hence we have an x* G (x',x r ) such that 

(x*) = 0. 

The unstable manifold W s (L(x)) is ^-symmetric to W u (L(x)). This means that 
a fragment of intersection of W s (L(x)) with S is parameterized by 

«4:[x',xn^M 2 , 
i4(x) = (tt x ™£ (x) ,-n Px w% (x)) . (7.12) 

Since w;^ (x*) = iu|; (x*) manifolds W"(L(a;)) and W s (L(x)) intersect at 

q* =g(x*,w u (x*)). 

Now we turn to proving transversality of the intersection at q* . By (7.10), around 
q* the manifold M(H(L(x))) is parameterized by x,y,p x - Therefore in the proof of 
transversality we restrict to these coordinates. Since Q is well defined, W u (L(x)) must 
transversally cross {y = 0}. By symmetry so does W s (L(x)). We therefore only need 
to prove that (x) and (x) intersect transversally in M 2 . 

Let x+ G (x*,x r ], 7 = l/(x+ -x*) and 

v = 7 ((x+, w"(x+)) - (x*, w"(x*))) G V+. 
By the mean value theorem 

wl (x+) - wl (x*) G n*# m -[DQ{B E )\ v. 

7 

By (7.9) this implies that 

7r a («£(x+)-«$(x*)) >0, 7T^K(x+)-^(x*))>0. (7.13) 
By mirror arguments, for x~ G [x , x*) 

ir x (w% (x~) - u£ (x*)) < 0, ^K(x-)-^(x*)) <0. (7.14) 
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From (7.13), (7.14) and (7.12) we see that w£ (x) and (x) intersect transversally 
at (x*) = wf, (x*) , which concludes our proof. □ 

Remark 7.9. From proof of Lemma 7.8 follows that we have the following esti- 
mate on the slope of the curves w^, (x) 

^pg(B E )v- - 

ir x DG(B E )V- ' 

By S-symmetry of W u (L(x)) and W s (L(x)) the slope of (x) is in —a. 
Once we verify (7.8) then by checking that a > we know that assumption (7.9) 
needs to hold. 

7.4. Proof of Theorem 4.1. In this section we write the computer assisted 
rigorous bounds, which we obtain using the method from sections 7.2, 7.3. As a result 
we obtain rigorous bounds for the position of fibers of W u (L(x)) and for transversal 
intersection of W u (L(x)) with W s (L{x)). By this we obtain the proof of Theorem 4.1. 

7.4.1. Bounds for Unstable Fibers. We start by writing out our changes of 
coordinates needed for application of Lemma 7.3 to the map (7.3) from section 7.2. 

We first choose the point qo = (x° ,0,0, p°) with x°,Py given in (4.1) and (5.5) 
respectively, i.e. 

x° = -0.9510055339445208, 
p y = -0.8368041796469730. 

We choose a matrix C as 

/ 0.197841 
_ -0.221884 
° - 1 

\ -0.255717 

We then choose four polynomials 

K Q (x) =0.1x - 0.0621591x 2 + 0.0375888x 3 - 0.0200645x 4 
K-l(x) =0.000533561x 2 - 0.00723085x 3 + 0.00827176x 4 
A' 2 (x) = - 0.0151949x 2 + 0.009304476x 3 - 0.00427633x 4 
K 3 (x) =0.0269670x 2 - 0.0275820a; 3 + 0.0203022x 4 

which define the nonlinear change of coordinates ip (see (7.2)). All of the above 
choices are dictated by (non-rigorous) numerical investigation. Above choice ensures 
that Ci(j(x, 0) + go gives a decent approximation of the position of the unstable fibers 
of q(x) for x E I for / given in (4.1). 

Now our computations start. We first compute the interval enclosure T such that 
r(q(x)) £ T for all x € I. The obtained result is 

T = [3.058882598,3.058883224]. 

Next we compute a set Bq such that (see Figure 7.2) 
1 p- 1 (C- 1 (q(x)-q ))(lB . 



TT Px pg(B E )v+ 

tt x DG(B e )V+ 



-0.197841 0.221884 \ 

-0.221884 0.773671 

1 -1 

0.255717 -1 1 
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Such set can be obtained using a technical Lemma 10.5 included in the Appendix. 
We thus obtain 

12,7.91575e- 12] \ 
12,7.91575e - 12] 
19,9.29424e - 19] 
08,4.50827e - 08] j 

Remark 7.10. Note that the set is flat along the third and stretched along the 
last coordinate. This is because we set up C and ip so that the third coordinate is 
associated with the section {y = 0} (on which lie q(x)) and that the last coordinate is 
associated with the direction of the curve q{x) = (x, 0, 0, k(x)). 

We now choose the size of our investigated set B in local coordinates and choose 
the parameters for our cones (see Figure 7.2). We take 

a = 2.56- 10~ 6 , 

and consider only one branch of the unstable manifold considering 

B= |J Q+(v)n{xe [x,x]} (7.15) 

veB 

with 

x=-l-10- n , x = 4.5-10~ 6 . 

The choice of x is dictated by the size of the fiber we later need to consider to 
prove intersections of stable/unstable manifolds. 

To compute a rigorous enclosure of [DF(B)] using (7.7), we subdivide B into 
N = 1200 parts Bi along the x coordinate 

B=\jB t . 

i=l 

Using Lemma 7.3 to obtain enclosures of F(Bi), combined with (7.7) we compute 
estimates for [DF(Bj)\. Combining the estimates [DF(Bi)} we obtain the following 
global estimate for [DF(B)] (the result is displayed with very rough accuracy, ensuring 
true enclosure in rounding) 



/ [-7.91575e- 
[-7.91575e - 
[-9.29424e - 

\ [-4.50827e- 



[DF(B)\ = 

/ [1465.6, 1466.5] [-0.353,0.369] [-0.285,0.283] [-0.300,0.333] \ 

[-0.361,0.360] [-0.360,0.361] [-0.290,0.277] [-0.319,0.304] 

[-0.138,0.140] [-0.139,0.139] [0.896,1.120] [0.458,0.700] 

\ [-0.201,0.202] [-0.202,0.202] [-0.171,0.149] [0.823,1.172] / 

Finally, using \DF(B)] and Lemmas 10.1, 10.2 we verify assumptions of Lemma 
6.3. We thus obtain rigorous bounds for the position of the fibers. The computation 
of the enclosure of the fibers took 18 minutes on a standard laptop. 

We plot the obtained bounds on fibers transported to the original coordinates 
of the system x,y,p x ,p y in Figures 7.3, 7.4, 7.5. On the plots we present rigorous 
enclosures of three fibers starting from q(x) with x on the edges and the middle of 
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interval / (with I chosen in (4.1)). This gives us an overview of the size of our fiber 
enclosures (left hand side of Figures 7.3, 7.4, 7.5). We can see that close to the set 
which contains {q(x) = (x, 0, 0, k(x))\x € /}, which is depicted in green, the estimates 
on the fibers is sharp (right hand plots in Figures 7.3, 7.4, 7.5). We can see that our 
three considered fiber enclosures are very close to each other, but are still separated, 
which is visible after closeup on the left hand side plot in Figure 7.5. 




Fig. 7.3. Projections of fiber enclosures onto x,y coordinates. 




-0.95100550 -0.95100545 -0.951005535 -0.951005530 



Fig. 7.4. Projections of fiber enclosures onto x,p x coordinates. 




Fig. 7.5. Projections of fiber enclosures onto x,p y coordinates. 

Remark 7.11. The range of obtained fibers is small. It is possible to reach 
somewhat further from q{x), but this significantly increases the time of computation, 
since further subdivision of the set is required. 

Remark 7.12. By using linearization only we have not been able to obtain ac- 
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curate enough enclosure of the fibers to handle the proof of transversal intersection 
of manifolds which follows in section 7.^.2. Thus the use of higher order change of 
variables seems to be needed. 

7.4.2. Bounds for Intersections of Manifolds. In this section we present 
rigorous-computer-assisted results in which we verify assumptions of Lemma 7.8 and 
thus conclude the proof of Theorem 4.1. For each x € I there are four points of 
intersection of W w (L{x)) and W s (L(x)) on {y = 0}. They can be seen on the left 
hand plot in Figure 7.6. We consider only the point which is furthermost to the right. 

We define the set Be — [x ( ,x r ] x B c with x',x r chosen as 

x' n = 4.461867506615821 • 10~ 6 , 

x ! =x m -10 -11 , x r =x m + 10- u . 

We verify that assumption (7.8) of Lemma 7.8 holds by computing Q(B l E ) and 
Q(B E ). We plot the obtained bounds in red on the right hand side plot of Figure 7.6. 
Next, using Remark 7.9 we compute 

a = [1.7695, 1.7725], (7.16) 

hence assumption (7.9) of Lemma 7.8 holds. Applying Lemma 7.8 concludes the proof 
of Theorem 4.1. 



0.4 | , , 1 1e-06 




-0.4 1 ' ' 1 -18-06 u u 

0.6 0.8 0.6500740 0.6500745 

X X 



Fig. 7.6. Left: Numerical sketch of W s (L(x)) (in red) and W u (L(x)) (in blue) intersected with 
{y = 0}- Right: we have proved that {ftx,p !C {W u (L(x")) n {y = 0})|x 6 /} consists of curves which 
pass through two red boxes. We have also proved that their slope is between [1.7695,1.7725], The 
blue/green line is a non-rigorous plot of the curves. 

We needed to subdivide Be into 600 parts to compute [DQ(Be)V + ] together with 
[DQ(Be)V~] with sufficient accuracy to obtain (7.16). Verification of assumptions of 
Lemma 7.8 took 24 minutes on a standard laptop. 

Remark 7.13. From a and by S- symmetry of manifolds W u (L(x)) andW 3 (L(x)), 
we obtain an estimate [58.8637°, 58.9439°] on the angle of intersection of the curves 
on the x,p x plane. 

Remark 7.14. At one go we obtain an estimate for a whole family of curves on 

Tr x , Px (W u (L(x)) (l{y = 0}) for x e /. 

In reality these curves are very close to each other (at furthest distance along x of 
about 2.65 ■ 10 _9 j. We plotted (using non-rigorous computations) two curves which 
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are furthest from each other on the right hand side plot of Figure 7. 6. One is in green 
and the other in blue. They are visible only after a large magnification, and on a 
paper printout will merge together. This means that our estimate on the position of 
the curves is somewhat rough in comparison to non-rigorous numerical simulation. 

8. Closing Remarks and Future Work. In this paper we have presented a 
method for proving existence of families of Lyapunov orbits in the planar restricted 
circular three body problem. The method gives explicit bounds on a curve of initial 
points, which can continue up to half the distance from L2 to the smaller primary in 
the Jupiter-Sun system. 

We also presented a method of proving transversal intersections of invariant man- 
ifolds associated with Lyapunov orbits. The method gives explicit bounds on where 
the intersection takes place. It has been applied to Lyapunov orbits with energy of 
the comet Oterma in the Jupiter-Sun system. 

In this paper we have focussed on detection of homoclinic intersections. Using 
identical tools one could also prove heteroclinic intersections of manifolds in the spirit 
of the work of Wilczak and Zgliczyhski [19, 20]. 

Due to the fact that the presented method gives explicit estimates on the position 
of investigated manifolds, it is our hope to later apply it to the study of diffusion. Here 
is an outline of future scheme that could be followed to prove diffusion. The family of 
Lyapunov orbits is normally hyperbolic, hence survives time periodic perturbations. 
In non- autonomous setting the system no longer preserves energy, which allows for 
diffusion between orbits of different energies. Such mechanism has been investigated 
in [G] for the planar restricted elliptic three body problem, for the system with special 
restriction on parameters. The discussed diffusion follows from the geometric method 
of Delshams, de la Llave and Seara [8, 9, 10] and requires computation of Mclnikov 
type integrals along homoclinic orbits of the PRC3BP. Since our method allows for 
precise and rigorous estimates for such orbits, it is our hope that such integrals could 
be computed using rigorous-computer assisted techniques. This combined with topo- 
logical methods [5, 7] for detection of normally hyperbolic manifolds could give first 
rigorous results for diffusion in the three body problem with real life parameters. 
From this perspective, the results of this paper are a first step in a larger scheme for 
investigation of real life systems. 

9. Acknowledgements. The author would like to thank Rafael de la Llave for 
discussions and remarks regarding implementation of the parameterization method. 
Special thanks go also to Daniel Wilczak for discussions on rigorous-computer-assisted 
computations using the CAPD library (http://capd.ii.uj.edu.pl). 

10. Appendix. 

10.1. Verification of Cone Conditions. Lemma 10.1. Let A be an interval 
matrix of the form 




where an = [o u , an] with a n > 0. If for any e <G e, \\e\\ < e holds 




(10.1) 



then for v = 
\\Av\\ > m \\v\\ 



(x, y) such that Q(v) = ax. 2 — ||y|| 2 > and any A £ A we have 
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Proof. For v = (x,y) satisfying Q(v) > 0, wc have ||x|| 2 + ||y|| 2 < ||x|| 2 (1 + a). 
Using (10.1) this gives the following estimate 

H^MI > «n ||x|| - e ||y|| > («n - (a/~o) ||x|| > ||x|| 2 + ||y|| 2 = m ||w|| . 

□ 

Lemma 10.2. Lei Q(i>) = <5(x,y) = ax 2 — ||y|| 2 , let Cq be a diagonal matrix 
such that v t Cqv = Q(v), and let A = [DF(Q + (v*))]. Assume that D = A t CqA is 
an interval matrix of the form 

Assume that dn = [d 1:L ,<iii] with d n > and that for some M > 0, for any sym- 
metric matrix B G B 

inf {A|A G spec (£?)} > -M. (10.2) 

If for any e G £ we /icroe ||e|| < e awe? d n — 2e > A/a, £/ien /or any v\, V2 £ U, vi =/= v-i 
such that Q (vi — v%) > 

Q - F(t> 2 )) > 0. 

Proof. By (6.9) Q (-F(vi) — F(v2)) = (vi — V2) T D (v\ — v 2 ) for some symmetric 
matrix DtD. 

For v = (x, y) such that Q( x > y) > and for any symmetric D E T> 

B 

we compute the following bounds 

v T Av = dux 2 + xe T y + y T ex + y T -By 

>d n X 2 -2e||y|||x|-M||y|| 2 

> (d n -2e)x 2 -M||y|| 2 

>MQ (x,y) 

> 0. 

□ 

Remark 10.3. Assumption (10.2) is easily verifiable from Gershgorin theorem. 

10.2. Bounds for the Images in Local Coordinates . Here we give a proof 
of Lemma 7.3. 

Proof. Inclusion (7.6) is equivalent to showing that for any r G T, and any Vi G U\ 
there exists an vi in U2 such that 

G(t, Vi ,v 2 )=0. (10.3) 
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Let us fix a t G T and v\ G U\ and use a notation G T _ vl {v 2 ) = G(t, Ui,^). Observe 
that [DG T>V1 (U 2 )] C A(U 2 ) and [G r>vi (v )] C [G(T, Ui, v Q )] . Since from (7.5) 

v - [DGr^iUa)]- 1 G T , Vl (v ) c N(T, v , Ui, U 2 ) C f/ 2 , 

by the interval Newton method (Theorem 2.1) there exists a unique v 2 = v 2 (r, ui) G 
U 2 , which satisfies (10.3). □ 

Remark 10.4. When applying Lemma 7.3, due to very strong hyperbolicity of 
the map $ T it pays off to use the mean value theorem. Taking U\ = v\ + B we can 
compute 



iV(T> ,[/i,[/ 2 ) 



A((7 2 )- 1 ||(T,[/ 1 ,«o) 



This is a better form since in (below we neglect arguments in order to keep the formula 
compact) 



_ x dG 
dv\ 



-DX 



1 ■ (Wr 1 



(10.4) 



the strong hyperbolic expansion cancels out. This is the main advantage of Lemma 
7.3. 

Here we give a technical lemma that can be used for computation of 

^(C -1 ((a, 0, (),«(&)) - q )) for x € I 

and q° = (x°, 0, 0,Py). Below, R can be any matrix close to Dip~ 1 (0)C~ 1 A. 
Lemma 10.5. Let a G R and J\ C R be from Lemma 3.1 and let 



A = 



( 1 \ 

10 

10 

\ a 1 J 



Let B be a set in M 4 , let R be a 4 x 4 matrix and let 

M := (i?B) i?] _1 (7 - a; , 0, 0, Ji - p°) . 

//' 

McB 

iften V^C" 1 ((a, 0, 0, k(x)) -q )) C i?B. 
Proof. By Lemma 3.1 

(s, 0, 0, k(x)) G 0, 0,p°) + (7 - x°, 0, 0, o(7 - x°) + J x - p°) 
= q° + A(L-x o 1 Q,0,J 1 -p o y ) 1 



(10.5) 



hence 



Let 



(s, 0, 0, k(x)) -q°eA(l- x°, 0, 0, Ji - pj) 



(10.6) 
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If we can show that for any q 6 [I — x° , 0, 0, J\ — p°) there exists apeB such that 

G q ( P )=0 (10.7) 

then 

V" 1 (C^Ag) = Rp, 

hence by (10.6) 

t/j^iC' 1 ((x,0,0,K(x))-q°)) c KB. 

To show (10.7) wc apply the interval Newton method (Theorem 2.1). Since V'(O) = 
we can compute 



/V(0,B) 



d 



G,(0) 



= - [A- 1 CL»^(i?B)i?] 1 (-q) 
C M, 

and by (10.5) combined with Theorem 2.1 obtain (10.7), and hence obtain our claim. 
□ 
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